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Abstract. Avalanches, or Avalanche-like, events are often observed in the dynamical behaviour of many 
complex systems which span from solar flaring to the Earth's crust dynamics and from traffic flows to 
financial markets. Self-organized criticality (SOC) is one of the most popular theories able to explain this 
intermittent charge/discharge behaviour. Despite a large amount of theoretical work, empirical tests for 
SOC are still in their infancy. In the present paper we address the common problem of revealing SOC from 
a simple time series without having much information about the underlying system. As a working example 
we use a modified version of the multifractal random walk originally proposed as a model for the stock 
market dynamics. The study reveals, despite the lack of the typical ingredients of SOC, an avalanche-like 
dynamics similar to that of many physical systems. While, on one hand, the results confirm the relevance 
of cascade models in representing turbulent-like phenomena, on the other, they also raise the question 
about the current state of reliability of SOC inference from time series analysis. 
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Econophysics 
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1 Introduction 

The theory of self-organized criticality (SOC) has been 
developed in the late eighty's by Bak, Tang and Wiesen- 
feld [1] , in order to explain the ubiquity of power laws in 
nature. The key concept of SOC is that complex systems 
"naturally" self-organize to a globally stationary intermit- 
tent state in which avalanche-like events are power law 
distributed. These features are similar to those found in 
physical systems at the critical point [2|. The prototyp- 
ical model of a system exhibiting SOC behaviour is the 
2D sandpile [T]. Here the cells of a grid are filled by ran- 
domly dropping grains of sand (external driving). When 
the gradient between two adjacent cells exceeds a certain 
threshold a redistribution of the sand occurs, leading to 
more instabilities and further redistributions. The char- 
acteristics of this system, indeed of all systems exhibiting 
SOC, is that the distribution of the avalanche sizes, their 
duration and the energy released, obey power laws. 

Remarkably, this kind of scale-free intermittent evolu- 
tion is similar to that observed in many physical and social 
systems. Examples include astrophysical and geophysical 
plasmas [Hll], earthquakes [5], evolutions of species [6], 
traffic dynamics [7], wars and the stock market [MTUl 
fTTl[T2] (a recent review on the subject can be found in 
Ref. [9]). Despite the theoretical interest, reliable tests to 
prove the presence of SOC in real systems are still in their 



infancy. Some attempts have been made in the contest of 
solar flaring [T^ , astrophysical P^^fTS] and laboratory (THl 
[TT] plasmas and the stock market [18]. These works, while 
leaving open the question of a SOC behaviour, clearly 
show that the evolution of these systems can be well de- 
scribed by an avalanche-like dynamics characterized by 
power laws in the avalanche size, duration and waiting 
time between them. The presence of correlation between 
laminar times, that is the time elapsed between two avalanches, 
in particular, has raised objections to the relevance of SOC 
in these contexts. In fact, due to the lack of memory in 
the random external driving commonly used in the simu- 
lations of conservative SOC systems, the probability dis- 
tribution function (PDF) of laminar times actually follows 
an exponential decay [19] . However, for non-conservative 
systems, power laws can still be observed in the presence 
of temporal correlations of the avalanches near the SOC 
state [201,1211 . Such temporal correlation could also be due 
to the intrinsic memory process (possibly chaotic) in the 
driver l22ll23l. 
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Motivated by recent observations of avalanche-like dy- 
namics in financial time series [18j , we investigate a possi- 
ble similar behaviour in the popular multifractal random 
walk (MRW) originally proposed in Ref. [24] as a paradigm 
for the stock market behaviour. This model, although not 
presenting the characteristic mechanisms of SOC, such as 
a threshold triggering for the avalanches, is able to repro- 
duce most of the stylized features of the stock market. 
Moreover, the MRW belongs to a family of cascade-like 
models widely used to reproduce the statistical features of 
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the velocity fluctuations in hydrodynamic turbulence and, 
therefore, the discussions outlined in the next sections go 
beyond their application to finance but can be extended 
to every complex system which displaying a turbulent-like 
dynamics. 

In the next section we introduce the asymmetric MRW 
proposed by Chen, Jayaprakash and Yuan |25j , which avalanche 
dynamics is investigate in detail in the rest of the paper. 
In Sec. [3] we introduce the method of analysis while the 
results are exposed in Sec. ID Discussions and conclusions 
are left for the last section. 



one-step volatility. The dynamics of the model, and there- 
fore the capability to reproduce the stylized facts of the 
stock market, is related to the dynamics of the variable 
St: small variations lead to an intermittent behaviour in 
r(t), similar to the one observed in the financial markets. 
Specifically, we can write St as 



St = So r 



(2) 



where Sq is related to the amplitude of the fluctuations 
while 7 to their intermittency. The term n{t), the core of 
the model, is a bounded random walk with increments 



2 The MRW Model: the CJY Version 



An{t) ^rjt+a W{t) - Pfj. 



(3) 



Here 



Recently, the study of the stock market, seen as a complex 
system, has attracted the attention of many physicists (for 
reviews see Refs. pSll^Tll^lT^ ). Its dynamical behaviour 
is characterized by "stylized facts" mainly concerning the 
logarithmic returns, r{t) = In [P{t)/P{t — 1)] (where P{t) 
is the stock price), and their absolute values, that can 
be regarded as a measure of the volatility, v{t) = |r(i)|. 
Such stylized facts have been used as a guide for validat- 
ing phenomenological models of stock price fluctuations. 
Among them, the appearance of "fat tails" in the PDF of 
the logarithmic returns, related to frequent large fluctu- 
ations in price, and the long time correlations present in 
the volatility, a phenomenon known as volatility cluster- 
ing, have been extensively investigated in the econophysics 
literature [29 |[30ll3T| [32| . The evidence of leptokurtic distri- 
butions in flnancial time series leads to immediate analo- 
gies with the longitudinal fluctuations of turbulent flows 
where a similar dynamics is observed, although differences 
have been pointed out as well . Motivated by the afore- 
mentioned evidence, cascade models, originally developed 
to reproduce the characteristic features of intermittency 
in hydrodynamic turbulence, have been applied, with suc- 
cess, to reproduce some stylized facts of the stock market 
dynamic^j. 

In this framework, one of the most popular models 
is the MRW originally proposed by Bacry, Delour and 
Muzy [24j. The version that we are going to use in the 
present investigation has been proposed by Chen, Jayaprakash 
and Yuan [25], referred to as CJY in the rest of the paper. 
In the CJY model the returns are expressed by 



<P(t) = Ki r]t-KN,+irit^ 



No 

-E 

1=1 



[K,+i-K,]r^t^„ (4) 



with rjt independent random variables, with average fj, 
which assume the values -1-1 with probability p and — 1 
with probability 1 — p. In our simulations fj = 2p — I < 0. 
This term, that alone can reproduce volatility clustering, 
has been found to be necessary in order to reproduce some 
scaling properties of the conditional fluctuations observed 
in financial data [25]. The second term in the increment, 
if'(t), has the ability to recover the long-time correlations 
of the market volatility and, from now on, we will refer to 
it as the multifractal increment Its strength coefficient, 
a, is related to the degree of intermittency of the time se- 
ries and, therefore, to the time scale of process. The kernel 
used in Eq.lgl) for the convolution is Ki = 1/^/i and we 
fix the memory steps to Nc ~ 1000 in the simulations. 
The last term of Eq. controls the drift rate in n{t), 
instead, with strength (3, and, therefore, adds more fiex- 
ibility to the model. A more detailed discussion on the 
present model with a — P — can be found in Ref. |41j . 
The parameter 7 is fixed to 1.05 in all the simulations and, 
following Ref. 25J, it is linked to a and p according to: 



ln(7)^ 



P 



1 



r 



(5) 



(6) 



r(t) ^ Stzt, 



(1) 



where zt is a Gaussian random variable with zero mean 
and unitary standard deviation, while St represents the 



^ In hydrodynamic turbulence "cascade" refers to the flow 
of energy from the largest scales, where it is injected, toward 
the smallest ones where it is finally dissipated. In the market 
contest, instead, it is assumed that there exists a flow of in- 
formation among the different temporal scales adopted by the 
traders. A further discussion on the subject is outside the scope 
of this work but the interested reader can refer to the seminal 
book of Frisch [33] for a general review of cascade models in 
turbulence and Refs. [34 
stock market. 



These particular choices are appropriate for a correct re- 
production of the observed statistical features of the data. 
Note also that, for the previous choice of the parameter 7, 
namely 1.05, the probability distribution of the random 
variable in the multifractal increment, Eq.((5]), becomes 
slightly asymmetric, p — 0.4756, in contrast with the sym- 
metry of the original MRW [24] . As previously mentioned, 
we impose reflecting boundaries for n{t) in order to pre- 
vent the realization of extremely large or small (unrealis- 
tic) fluctuations in the simulation of the market activity. 



for applications to the 



^ Eq. (13} , is related to the original formulation of the MRW in 
Ref. [24] where the logarithmic variance is expressed in terms of 
a convolution between a memory kernel and a random process 
as reported in Ref. [39ll40] . 
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Fig. 1. (a) Part of the time series of logarithmic returns for 
daily data of the S&P500 index, used to calculate the PDF in 
(c). The all set spans from 3/1/1950 to 18/7/2003. (b) Time 
series generated from the model of Eqs. ((2)| and (|3]) with pa- 
rameters 7 = 1.05, ao = 0.1, P — 1.3 and So ~ 1. The two 
time series shown in the plot have been standardized accord- 
ing to r{t) —> {r{t) — {r))/a{r), where (. . .) and a represent, 
respectively, the average and the standard deviation over the 
period in examination, (c) Comparison between the PDFs in 
(a) and (b). The PDF generated with the parameters 7 = 1.05, 
ao = 0.2, /3 — 4.0 and (So = 1 is also shown. 



namely < n{t) < Umax, with Umax = ln(30)/ 111(7). By 
doing that, 5t in Eq.((T]) is bounded between 5q and 30 (5o. 

In Fig. [1] we compare the time series generated with 
the model of Eqs. dH) and Fig. □ (b), and the time 
series of daily returns for the S&P500 inde^ll, Fig. [1] (a). 
The parameters used in the simulation, and reported in 
the caption of the figure, have been chosen in order to 
match the properties of the financial data set, as under- 
lined by the similarities in the PDFs of the two processes. 
Fig. [1] (c). Note, however, some discrepancies in the tails 
of the PDFs. A possible explanation is that the dynam- 
ics of the extreme events differs from the dynamics of the 
bulk of the distribution (generated by the CYJ model, in 
this case) and they could be interpreted as "outliers" [42l 
I40j . However, finite size effect in the relatively short time 
series of the S&P500 should be also considered. 



3 Wavelet Transform Filtering and Analysis 
Method 

As mentioned in the introduction, many complex systems 
show an intermittent activity: quiescent periods are sud- 
denly interrupted by bursts of activity. This kind of non- 
stationary dynamics is often related to multi-scale phe- 
nomena and most of the time standard analysis techniques 
can fail to reveal some important events that are localized 
in time or scale [33]. This is, for example, when using 
elementary filters: along with the noise background also 
meaningful information can be filtered out [18j . 



3 The data have been collected from 3/1/1950 to 18/7/2003 
for a total of 13468 samples. 



In order to overcome these problems, wavelet based 
techniques are becoming more and more popular in com- 
plex systems applications [43j . This approach enables one 
to decompose the signal in terms of scale and time units 
and so to separate its coherent parts (or "avalanches") 
- that is, the bursty periods related to the tails of the 
PDF - from the noise-like background, thus enabling an 
independent study of the intermittent and the quiescent 
intervals plfTelHtlliSlfTSj . 

The idea behind the wavelet transform is similar to 
that of windowed Fourier analysis and it can be shown 
that the scale parameter is indeed inversely proportional 
to the classic Fourier frequency. The main difference be- 
tween the two techniques lies in the resolution in the time- 
frequency domain. In Fourier analysis the resolution is 
scale independent, leading to aliasing of high and low fre- 
quency components that do not fall into the frequency 
range of the window. However, in the wavelet decompo- 
sition the resolution changes according to the scale (i.e. 
frequency). At smaller scales the temporal resolution in- 
creases at the expense of frequency localization, while for 
large scales we have the opposite. For this reason the 
wavelet transform can be considered a sort of mathemat- 
ical "microscope". While the Fourier analysis is still an 
appropriate method for the study of harmonic signals, 
where the information is equally distributed, the wavelet 
approach becomes fundamental when the signal is inter- 
mittent and the information localized. 

For a time series analysis it is often preferable to use a 
discrete wavelet transform (DWT) . The DWT can be seen 
as a appropriate sub-sampling of the continuous wavelet 
transform (CWT) by using dyadic scales. That is, one 
chooses A = 2-', for j — 0, L — 1, where L is the number 
of scales involved, and the temporal coefficients are sep- 
arated by multiples of A for each dyadic scale, t = n2^ , 
with n being the index of the coefficient at the jth scale. 
The DWT coefficients, Wj,n, can then be expressed as 

W^,,n = (/, ^,,n> = 2-^/2 J f{u)^j{2~^u - n)du, (7) 

where ipj^n is the discretely scaled and shifted version of 
the mother wavelet. The wavelet coefficients are a mea- 
sure of the correlation between the original signal, /(<), 
and the mother wavelet, ipi^) scale j and time n. For 
the DWT, if the set of the mother wavelet and its trans- 
lated and scaled copies form a complete orthonormal basis 
for all functions having a finite squared modulus, then the 
energy of the starting signal is conserved in the wavelet 
coefficients. This property is, of course, extremely impor- 
tant when analyzing physical time series |45| . Following 
Ref . [18^ , we use the Daubechies-4 as mother wavelet [46] 
for the analysis presented in the next sections. Tests per- 
formed with different sets do not lead to any qualitative 
difference in the results. A more comprehensive discus- 
sions on the general properties of wavelets and their ap- 
plications are given in Refs. [151133] • 

The importance of the wavelet transform in the study 
of turbulent-like signals lies in the fact that the large 
amplitude wavelet coefficients are related to the extreme 
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Fig. 2. PDFs of the original time series, Fig. [3] (a), obtained 
with the model of Eq. (O (7 = 1.05, ao = 0.1, /3 = 1.3 and 
So = 1) and its filtered version. Fig. |3](b). A Gaussian is also 
plotted for visual comparison (dashed line) . The Daubechies-4 
wavelet used for the analysis is shown in the inset. 



events in the tails of the PDF, while the laminar or qui- 
escent periods are related to the ones with smaller ampli- 
tude |45| . In this way it is possible to define a criterion 
whereby one can filter the time series of the coefficients 
depending on the specific needs. In our case we adopt the 
method used in Refs. [iSllTlilTB] and originally proposed 
by Katul et al. [l^. In this method wavelet coefficients 
that exceed a fixed threshold are set to zero, according to 



if 



otherwise. 



wl^ < C 



(8) 



here (...)„ denotes the average over the time parameters 
at a certain scale and C is the threshold coefficient. Once 
we have filtered the wavelet coefficients Wj^n we perform 
an inverse wavelet transform, obtaining a smoothed version 
of the original time series. 

The residuals of the original time series with the fil- 
tered one correspond to the bursty periods, or avalanches, 
which we aim to study. An example of the filtering tech- 
nique in terms of PDFs is given in Fig. [2l 

Once we have isolated the noise part from our signal 
series we are able to perform a reliable statistical analysis 
on the avalanches of the residual time series. In particular, 
we define the avalanches as the set of events during which 
the volatility of the residual time series, Vres{t) = |rres(i)|, 
is constantly above a positive small threshold, e « 0. It 
is also possible to find an optimal value for the choice 
of the threshold parameter C in a way that the filtered 
time series is, as close as possible, uncorrelated Gaussian 
noise. From previous studies, this parameter is found to be 
C ~ 1 [18]. In any case, the resulting statistical analysis 
is qualitatively unchanged as long as < C < 4 [45lfT8]. 
A graphic example of the procedure for extracting the 
avalanches is illustrated in Fig. [3l 

In analogy with the dissipated energy in a turbulent 
flow, we define the size of an avalanche, E, as the inte- 
grated value of the squared volatility over each coherent 
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Fig. 3. (a) Window of the time series obtained with the model 
of Eq. © (7 = 1.05, ao = 0.1, /? = 1.3 and So = 1). (b) Filtered 
version of the time series in (a) obtained with the Daubechies-4 
wavelet. In order to properly identify the avalanches we create 
a residual time series by subtracting the filtered time series 
from the original one. This noise removal technique presents 
advantages to standard threshold methods when applied to 
multiscale systems [18j . (c) Absolute value of the residuals, or 
residual volatility, used to extract the avalanches, as explained 
in the text. Two examples of avalanches, among the many pre- 
set, are underlined by the dashed lines in (c). The horizontal 
line slightly above zero represents the small threshold e. 



event of the residual time series. The duration, D, is de- 
fined as the interval of time between the beginning and 
the end of a coherent event, while the laminar time, L, 
is the time elapsing between the end of an event and the 
beginning of the next one. 



4 Time Series Analysis 

In the previous section we have seen how the wavelet 
multi-scale filtering technique is an excellent tool to re- 
move uncorrelated Gaussian noise from an input signal. In 
particular, this filtering method becomes relevant when- 
ever the examined time series presents an intermittent be- 
haviour, that is an irregular switching, between periods 
characterized by large fluctuations and noise-like ones. By 
using this technique, the avalanches, which characterize 
an emergent behaviour in the dynamics, are highlighted at 
the expense of the uninteresting background. In this way 
we can make a proper statistical analysis of the quantities 
that characterize these coherent events, Fig. |31 

The statistical study of the avalanches identified with 
the wavelet technique is of great interest, not only because 
this would further test the capability of this model to re- 
produce the styfized facts of the stock market, but it could 
also shed some light on the relevance of this test in distin- 
guishing between SOC and non-SOC processes in a time 
series analysis. 

The analysis is carried out by studying how the statis- 
tics of the avalanches change as we tune the parameters of 
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the CJY model. The time series generated with this algo- 
rithm have a length N ^ 5 ■ 10^ {N = 2^^). Moreover, we 
investigate the specific relevance of each term composing 
the increments of the variable n{t), Eq.([3]). In order to do 
this, we independently analyze the time series generated 
with different expressions for the increments, An{t): for 
each simulation we consider a random walk for n{t) with 
boundaries rimin = and Umax = ln(30)/ ln(7). 

We also fix C = 1 as the threshold coefficient for the 
wavelet analysis. This particular value of C is close to the 
optimization value in the de-noise procedure. However a 
different choice of this parameter would not change the 
qualitative results of our analysis [H] . 

4.1 The Role of Multifractal Increments 

We first investigate the avalanche dynamics generated by 
the multifractal increments, that is the part of the CJY 
model that is related to the original formulation of the 
MRW. In this case Eq.^ reads as 

An{t)^aW(t), (9) 

where 'F{t) is given by Eq.ljl]) and the strength coefficient 
is expressed by a = ao/ ln(7) with 7 — 1.05. 

The avalanche analysis, resulting from the wavelet fil- 
tering, for the size E, duration D, and laminar times L, 
is carried out for different values of uq. This parameter, 
and a as a consequence, is related to the degree of in- 
termittency of the time series and, therefore, to the time 
scale of the process. For ao < 0.05, the dynamics of r{t) 
is dominated by noise and its PDF is a Gaussian. In the 
stock market contest, as well as in turbulence, this corre- 
sponds to observe the fluctuations at large scales. As we 
move this parameter toward larger values, the time series 
of r{t) becomes more and more intermittent, giving rise 
to the large fluctuations which characterize the broad tails 
of the PDFs of turbulent phenomena at small scales. The 
results are shown in Figs. [H [5] and [6] 

From these plots it is possible to observe how the 
purely multifractal model is not able to reproduce the 
power law behaviour observed in real data sets. In par- 
ticular, the PDFs for E and D of the avalanches decay 
exponentially, independently on the value of ao. This is an 
indication of the randomness behind the avalanche gener- 
ation process. The distribution of laminar times in Fig. [6l 
instead, show a Poisson-like shape for small values of this 
parameter, ao = 0.05, while they start slowly to converge 
toward a power law shape, P{L) ~ L~'^, for ao > 0.2. The 
resulting exponent, v ~ 2.3, is similar to the one found in 
the empirical studies jT8]. 

For completeness, we study also the avalanche dynam- 
ics by using the standard implementation of the MRW [39J . 
In this case 5t = (5oe"^*^ and An{t) is given by Eq. ^ 
where, this time, 77^ is a Gaussian random variable with 
zero mean and unitary standard deviation. The shapes of 
the PDFs of E and D resulting from the analysis, shown 
in the Figs. 2] and [5] for a = 0.25, display a similar a de- 
pendence as the CYJ version, as expected. However, the 
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Fig. 4. PDFs of the avalanche sizes, E, as a function of the 
parameter qq- For each value of this parameter the time series 
of E is noise-like and the correspondent distribution displays 
an exponential decay. The strength ao controls the decay rate 
of the distribution. The PDF for the MRW (X) with a = 0.25 
is reported as well. 
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Fig. 5. PDFs of the duration, D, for the purely multifractal 
model. No power law behaviour is observed and the same con- 
siderations as for the size, E, hold in this case. The PDF for 
the MRW (X) with a = 0.25 displays a similar behaviour. 

PDF of laminar times, at least at small temporal scales, 
FiglHl display a clearer power low behaviour compared to 
the CYJ model. 



4.2 The Complete CJY Model 

We turn now our attention to the complete CYJ model of 
Eq.(I3]). In this case, it has been shown [25] that a proper 
tuning of the parameters can reproduce most of the styl- 
ized features of the stock market. A particular good agree- 
ment between the model and the empirical data has been 
found by fixing 7 = 1.05 and Sq = 1.0, for the two couples 
of parameters (ao = 0.1,(3 = 1.3) and (ao = 0.2, /? = 4.0) 
as strength parameters for the memory term and the drift 
respectively. It is, therefore, of particular interest to ex- 
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L 



Fig. 6. Distribution of laminar times, L, in the multifractal 
model. A slow convergence from an exponential decay to a 
power law one is observed for ao > 0.2. The exponent of the 
power law, in this case, is u ^ 2.3. In this case the PDF for 
the MRW (X), with a = 0.25, displays a clearer power law 
behaviour compared to the CYJ version. 
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Fig. 7. PDFs of the size, E, in the case (ao = 0.1, P = 1.3) and 
(ao = 0.2, /3 = 4.0). A similar power law shape, with exponent 
V ~ 2.6, is observed for the two distributions. A dashed line is 
also plotted for visual comparison. 

plore the avalanche behaviour using these parameters. The 
results of the analysis are shown in Figs. [3 [8] and [9] for E, 
D and L. 

In this case, for some order of magnitude, we find a 
power law behaviour for the quantities under consider- 
ation. This is in qualitative agreement with the results 
found for the stock market. In particular, the exponents of 
the power laws seem to be close to the ones found for the 
analysis of the tick-by-tick Nasdaq E-mini Futures jl8. . 
Note that a scale-free avalanche dynamics has also been 
observed in other reduced models of turbulence, the shell 
models [13], via a simple threshold technique. 

This result can have important consequences regard- 
ing the possible identification of SOC in the stock market, 
and other complex systems in general, through a time se- 
ries analysis. In fact, we have shown that an avalanche- like 




Fig. 8. PDFs for the avalanches duration, D, for the same 
parameters of Fig. [T] In this case a power law with exponent 
v 4.7, is observed. 




L 



Fig. 9. PDFs for the laminar times between avalanches, L. 
The parameters used are the same as in Fig. [T] A power law 
with exponent u ~ 1.8 and a cut-off at L ~ lO'^ is reported. 

behaviour can also be observed in models, such as the one 
presented in this work, in which the characteristic ingre- 
dients of SOC, such as threshold dynamics, are actually 
missing. This, of course, does not rule out the possibility of 
SOC but, nevertheless, more relevant and discriminating 
tests become necessary. 

We further investigate our model by studying how the 
distribution of the E, D and L change with the drift 
strength /3. In fact, different markets could have differ- 
ent power law exponents, no universality has been found 
until now, and we want to test the elasticity of the CYJ 
with respect to this parameter. The PDFs for E, D and 
L as functions of /3 are reported, respectively, in Figs. fTOl 
[Illandini 

The analysis shows how the parameter /3 plays an im- 
portant role in the dynamics of the model. In fact, al- 
though the shape of the PDFs for E and D are robust 
against variations of the exponent changes with this 
parameter. Higher values of P imply a larger value for 
ly. The same arguments do not hold for the statistics of 
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Fig. 10. Dependence on the parameter /3, with uq = 0.1, for 
the PDF of the size, i?, of the avalanche. The drift strength 
controls the slope of the power law, which appear to go to 
saturation for the higher values of /3. 
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L 

Fig. 12. The exponent of the power law of the laminar times 
between the avalanches, L, do not seem to be sensitive to 
changes in the parameter /3, which appear to be relevant just 
in changing the cut-off of the distribution. 
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Fig. 11. Dependence on the parameter /3, with ao = 0.1, for 
the PDF for the duration D of the avalanches. Also in this case, 
the drift strength controls the slope of the PDF and saturates 
for large /3. 

the laminar times, L. In this case the resulting distribu- 
tion is pretty much independent of /3 and could constitute 
a limit in the present formulation of the model. A sep- 
arate discussion is reserved for /3 = 0. In this case the 
filtering procedure, with C = 1 as optimal value, is not 
able to remove all the wavelet coefficients related to the 
large fluctuations. As a consequence, the excess of kur- 
tosis, Ke = {r'^) / {r'^)'^ — 3, of the filtered time series, al- 
though still small in absolute value, K^. ~ 0.2, becomes 
more than one order of magnitude larger that in the pre- 
vious analysis . This means, to some extent, that there is 
not enough Gaussian noise to be filtered out in the time 
series! Moreover, the shapes of the PDFs related to the 
avalanche dynamics, show a behaviour that is systemat- 
ically different from the one observed once the drift is 
included. Figs. [lOl H and [H 

As a final note, we consider, for (3 — 1.3 fixed, how our 
analysis would change in the absence of the multifractal 
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Fig. 13. PDFs for E (Top), D (Middle) and L (Bottom) for 
(3 = 1.3 and ao ~ 0.1, with multifractality, and with qq = 0, 
without the multifractal term. 

increment, a^{t), in Eq.(I3]). In doing so, we compare the 
case with = 0.1 and ao = 0. The results are shown in 
Fig.Ill 

While the lack of the multifractal/memory term does 
not alter the distribution of laminar times between the 
clusters of volatility, it does increase the steepness of the 
distribution of E and L. This result is actually expected 
since this term builds the correlations inside periods of 
high volatility of the market time series, that is the so- 
called volatility clustering. By removing it we explicitly 
cut off part of the correlations inside the model, resulting 
in shorter avalanches. 



5 Discussion and conclusion 

In the present work we investigated a possible avalanche- 
like dynamics in an extended version of the popular multi- 
fractal random walk, the CJY, proposed as a paradigm for 
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M. Bartolozzi: Scale-free avalanches in the multifractal random walk. 



the stock market dynamics. We have been able to identify 
avalanche-like events in the fluctuations generated by this 
model. Subsequently, the statistical properties of these 
events have been estimated. The identification of these 
clusters goes through an intermediate passage where we 
use a wavelet filtering technique in order to suppress the 
contribution of the noise background and, therefore, en- 
hance the precision of our measures. 

The results show that, for a broad range of the param- 
eters, the distribution of size, duration and laminar time 
between avalanches follow a power law distribution. A very 
similar behaviour has been found in empirical studies on 
financial time series [18]. Therefore, we confirm the rele- 
vance of the cascade models, and in particular the CYJ 
version of the MRW, in modelling financial market. Our 
results also extend beyond the financial environment since 
this framework is quite general for describing dissipative, 
intermittent, systems, such as solar fiaring or MHD turbu- 
lence, where an avalanche-like dynamics has been observed 
as wefl pllMlfTSlfTBlfTT]. 

Equally important, our results stress how models lack- 
ing of typical SOC mechanisms can, nevertheless, mani- 
fest an avalanche- like behaviour. In fact, the recognition of 
SOC "patterns" could be an artifact of the identification 
method itself. In this regard, it is worth mentioning that 
simple threshold techniques, when applied to time series 
generated by SOC-free processes, can lead to the iden- 
tification of avalanches with size distributions (according 
to some specific definition) that are power laws. For ex- 
ample, it is well known that the "time for first return to 
the origin" of a random walk is power law distributed [40j 
despite the lack of any correlation in the time series. A 
more sophisticated example is reported in Ref. [l^ where 
it has been shown that the behaviour of a self-affine frac- 
tional Brownian motion, when analyzed with a moving 
average technique, can mimic the avalanche dynamics of 
the Dhar-Ramaswamy sandpile. These results illustrate 
the technical ambiguity in the identification of SOC from 
time series. As long as we do not have any a priori infor- 
mation about the underlying dynamics of the system, it 
is very hard to tell if the avalanches that we observe are 
a result of a genuine SOC dynamics or any other diffusive 
process. 

In conclusion, SOC has been claimed, perhaps too loosely, 
to play a role in many complex systems, but there is no 
method that is reliable enough to test its presence from the 
analysis of noisy time series. This is a very relevant issue 
since, in practical situations, all the available information 
regarding a system is encoded in its time series. Therefore, 
an extension of this theoretical framework, which would 
enable the present gap with empirical analysis to be filled, 
is of great practical importance and would probably settle 
many speculations on the subject. 
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